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ABSTRACT 

We present a set of hydro-dynamical numerical simulations of the Antennae galaxies 
in order to understand the origin of the central overlap starburst. Our dynamical 
model provides a good match to the observed nuclear and overlap star formation, 
especially when using a range of rather inefficient stellar feedback efficiencies (0.01 < 
9EoS ^ 0.1). In this case a simple conversion of local star formation to molecular 
hydrogen surface density motivated by observations accounts well for the observed 
distribution of CO. Using radiative transfer post-processing we model synthetic far- 
infrared spectral energy distributions (SEDs) and two-dimensional emission maps for 
direct comparison with Herschel-PAGS observations. For a gas-to-dust ratio of 62:1 
and the best matching range of stellar feedback efficiencies the synthetic far-infrared 
SEDs of the central star forming region peak at values of ~ 65 — 81 Jy at 99 — 116^m, 
similar to a three-component modified black body fit to infrared observations. Also the 
spatial distribution of the far-infrared emission at 70/im, 100/im, and 160/im compares 
well with the observations: > 50% (> 35%) of the emission in each band is concentrated 
in the overlap region while only < 30% (< 15%) is distributed to the combined emission 
from the two galactic nuclei in the simulations (observations). As a proof of principle we 
show that parameter variations in the feedback model result in unambiguous changes 
both in the global and in the spatially resolved observable far-infrared properties of 
Antennae galaxy models. Our results strengthen the importance of direct, spatially 
resolved comparative studies of matched galaxy merger simulations as a valuable tool 
to constrain the fundamental star formation and feedback physics. 

Key words: methods: numerical - galaxies: evolution - galaxies: individual: NGC 
4038/4039 - galaxies: interactions - galaxies: starburst. 



1 INTRODUCTION 

Being one of the closest major mergers, NGC 4038/39 (Arp 
244; also known as "the Antennae") is an ideal laboratory 
for studying the details of interaction-induced star formation 
in the local Universe. Especially the overlap region between 
the two interacting disks is surprisingly rich in molecular 



gas (Wilson et al. 2000 Gao et al. 2001 Zhu et al. 20031 



and is subject to strong large-scale shocks driven by the 
galaxy collision ( Haas et al.|2 005 Hcrr era et al.|201~T 2012 



* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and 
with important participation from NASA, 
f E-mail: skarl@ast.cam.ac.uk 



|Wei et aL||2012| |. The spatial distribution of the molecular 
gas as traced by CO emission correlates well with the large 
population of young massive star clusters in the Antennae 
(|Wilson et al.|2000| [Zhang et al.|200T||Whitmore et al.|2010| 
Ueda et al.|2012[ ). 



Consistently, early mid-infrared ISO observations 
showed that a significant amount of the total star forma- 
tion of the Antennae takes place in the overlap region ( |Vi- 
groux et al. 1996. Mirabel et al. 1998). These results were 



confirmed by Spitzer mid-infrared observations ( |Wang et al 
|2004b|> and recent Herschel-PACS far-infrared (FIR) data 



I Klaas et al. 20101. By combining different bands in the 



FIR, Klaas et al. (20101 determined local star formation 



rates (SFRs) in small, localized regions of the Antennae. 
These observations showed that the Antennae galaxies har- 
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bour a major off-nuclear starburst in the overlap region and 
an arc-like star forming area in the northern galaxy. This 
makes the Antennae an exceptional system, together with 
only a handful of other mergers showing similar levels of 
off-nuclear interaction-induced star formation, e.g. Arp 140 



[Cullen et al. 20071, II Zw 096 (Inami et al 



6090 (jDinshaw et al.|p99} |Wang et al.||2004a[ ), NGC 6240 



2010), NGC 



( Tacconi et aL|l 999; Eng el et al.|2010[ ), and NGC 2442 ( |Pan 
coast et al.|2010 \. This behaviour contrasts with the general 
observed trend in a sample of 35 pre-merger galaxy pairs 



observed by Smith et al. ( 2007 1 and with results from previ- 
ous numerical simulations (e.g. Barnes & Hcrnquist 19961. 
In the Antennae, the total SFR in the overlap region even 
exceeds the combined SFR of the two galactic nuclei by a 
factor of ~ 3 (Klaas et al. 20101. While most of the star 



forming regions across the Antennae show a modest level 
of star formation, a few areas in the overlap region and the 
western loop are sites of intense localized bursts with specific 
star formation rates similar to those of heavily star-bursting 
ULIRGs ( [Wang et al.|2004b| [Zhang et al.|2010" |. 



Until recently, numerical models were not able to re- 
produce the extraordinary star formation properties of the 
Antennae. They failed to produce a sufficiently gas-rich over- 
lap region and underestimated the level of total star for- 
mation in the Antennae, which is measured at a rate of 



~ 5 


- 22M yr" 1 (e.g. 


Stanford et al. 


1990| |Zhang et al. 


2001 


|Knierman et al.|2003; Brandl et al.|2009; Klaas et al. 


2010 


). However, in the last few years there have been a num- 



ber of promising numerical studies of the Antennae, e.g. 
about the nature of the observed star formation and the 
current interaction stage ( Renaud et al.|[2008| |Karl et al.| 



2010 Teyssier et al. 20101, the interpretation of the ob- 



served cluster age distribution (Karl et al. 20111, and the 
generation and evolution of magnetic fields in this proto- 
typical galaxy merger ( Kotarba et al.|[2010 l. The nature of 
the overlap starburst has been interpreted in two ways: (1) 
the gas-rich overlap region might originate from the two pro- 
genitor disks being observed in projection on the plane-of- 
the-sky (Barnes 19881. The high level of star formation in 
a spatially extended pattern may result from efficient frag- 



mentation and subsequent star formation (Teyssier et al 



2010); (2) the actively star forming overlap region may be a 
natural consequence of the special interaction phase of the 
Antennae, shortly before the two galaxies merge, at a time 
when both disks already overlap, interact and large-scale 



shocks induce localised star formation (Zhang et al. 2010 
Karl et al.|2010| |Whitmore et al.||2010[ ). 



In this paper we re-visit the role of the physical prop- 
erties in the star forming interstellar medium (ISM) for the 
merger-induced star formation in the Antennae. In particu- 
lar, we aim to constrain the net-effect of the applied thermal 
stellar feedback on the star formation properties at the time 
of best match in a number of simulations employing brack- 
eting cases of the parameter of the effective equation of state 
<?EoS in our adopted stellar feedback algorithm. We find that 
the star formation histories in the different simulations are 
complex and there is a large scatter in the global observable 
star formation properties. Additional information provided 
by a spatially resolved analysis of these properties therefore 
seems to be needed to break degeneracies in simulations. To 
this end, we post-process the simulation results of the An- 
tennae simulations, for the first time, using an efficient 3D 
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Figure 1. Observed integral spectral energy distribution of the 
Antennae, from ISOCAM 15 /im, Spitzer-MIPS 24 fim, Herschel- 
PACS 70, 100, and 160 /im, Herschel-SPIRE 250, 350, and 500 Atm 
and IRAM-MAMBO 1200 fim photometry. Black triangles are the 
flux values originating from the standard calibration of the respec- 
tive instrument. Red squares are colour corrected fluxes (colour 
correction only applied for the PACS and SPIRE photometry) 
according to the thermal components (dotted lines) fitted to re- 
sult in a composite model SED (solid black line). Details of the 
modified black body fits are given in the figure legend as well as 
the 850 fim flux prediction used for dust mass determination and 
the total 10 - 1000 fim flux. 



Monte Carlo radiative transfer method to compute the con- 
tinuum spectral energy distributions (SEDs) and to produce 
synthetic maps in different bands in the far-infrared, which 
are then compared to detailed, spatially resolved far-infrared 
observations of star formation sites in the Antennae ( |Klaas| 
et al.|20To see also Section 2.3 1. 

We describe the merger simulations as well as details of 
the radiative transfer post-processing procedure and the FIR 
observations in Section [2] In Section [3] we discuss the star 
formation properties in the simulations subject to varying 
efficiencies in the stellar feedback model. We compare to 
the FIR observations in Section [4] Finally, we conclude and 
discuss in Section [5] 



2 METHODOLOGY 

2.1 A-body+SPH simulations 

The numerical methods and simulations used in this study 
have, to a major part, been presented in previous papers 
( Karl et al.|2010||201l[ see also |Johansson et al.|2009| ). Here, 
we briefly summarize the most important facts and describe 
the details of the new simulation runs. All simulations were 
performed using the parallel Trec/SPH code Gadget-3 (see 



Springcl 2005 ) including cooling, star formation and thermal 
stellar feedback as detailed below. 

Following Kat z et al.| ( | 1996] ) , radiative cooling is com- 
puted for an optically thin plasma composed of primordial 
hydrogen and helium, assuming ionization equilibrium in 
the presence of a spatially uniform and time-independent 



local UV background (Haardt & Madau 1996). Star forma- 
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Figure 2. Projected gas surface density and star formation in the central ±18 kpc of the galactic disks for the three Antennae simulations 
(K10-Q1, K10-Q0.1 and K10-Q0.01) with different values for the parameter for the softened equation of state, qeoS = 1-0 (left), <7EoS = 0-1 
(middle), and QeoS = 0-01 (right). Top panels: Projected gas surface density overlaid with contours of the star formation rate surface 
density. The contours correspond to levels of log(Sgpri /[Mq yr _1 kpc -2 ]) = —6, —5,..., —3, —2, respectively. Bottom panels: Gas 
surface density contours and stellar particles formed within the last 100 Myr, colour-coded by their age (see legend). 



tion and feedback are treated using the effective multi-phase 



model by Springel & Hernquist (20031, where stars are al- 
lowed to form on a characteristic timescale t+ in regions 
where gas densities exceed a critical threshold density of 
"-crit = 0.128cm -3 . In these star forming regions, the in- 
terstellar medium (ISM) is assumed to develop an effective 
two-phase fluid, where cold clouds are embedded in a hot 
ambient medium in pressure equilibrium. This is achieved 
by instantaneously returning mass and energy from mas- 
sive short-lived stars to the ISM upon formation of stellar 
particles; we chose a mass fraction of f3 = 0.1 in massive 
stars (M+ > 8 Mq) and 10 51 ergs -1 of energy released per 
supernova, according to a Salpeter-type stellar initial mass 
function ( Salpeter|1955 1. 



As a result of the two-phase formulation, the "effec- 
tive" equation of state (EoS) in star forming regions is quite 
"stiff" , providing strong and steeply rising pressure support 



to the gas with increasing density (see Figure 1 in Springel 
fc Hernquist|2003[ ). A convenient way of parametrizing un- 
certainties in the complex sub-grid physics is to control the 
efficiency of the thermal feedback via a further dimensionless 



parameter, qeoS ( Springel et al.|2005 l, which linearly inter- 
polates between a full, "stiff" feedback (qeoS = 1) and a 
softer, isothermal equation of state with T = 10 4 K (qsos = 
0). The resulting equation for the pressure in star forming 
regions then takes the form P — qeoS Pea + (1 — <7Eos) -Rso, 
where P cfi = (7 - 1) pu cft and P iso = (7 - l)p« iso . In 
this case, 7 = 5/3 is the adiabatic index of an ideal gas, 
Ueg = zitcoid + (1 — x)uixot, the effective specific internal 
energy weighting the hot and cold phases by x, the mass 
fraction in cold clouds, and Ui so is the specific internal en- 
ergy corresponding to T = 10 4 K. In this study, we test 
our Antennae model with various choices for qeoS against 
changes in the quality of the match to Herschel-PACS ob- 
servations in the FIR. 

Motivated by comparable total magnitudes in the B- 



and R-band for both NGC 4038 and NGC 403 9 ( [Lauberts 
|fc Valentijn||1989| |de Vaucouleurs et al.||1991[ ), we model 
the Antennae system as a 1:1 merger with a total mass of 
Mtot = 5.52 x 10 11 Mq for each galaxy. The initial disk mod- 



els are set up following the method described in |Springel| 
et all (|2005h. They all consist of a iHernquistl (|1990|) cold 
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Table 1. Nomenclature and discerning characteristics for the numerical simulations. 



Mode[°] 






-^disk.gas 


Afbulge 


9EoS 


k 


RT eff. resolution 


K10-Q0 


55.2 


3.3 


0.83 


1.4 





20% 


2048 3 


K10-Q0.01 


55.2 


3.3 


0.83 


1.4 


0.01 


20% 


1024 3 


K10-Q0.1 


55.2 


3.3 


0.83 


1.4 


0.1 


20% 


2048 3 


K10-Q0.5 


55.2 


3.3 


0.83 


1.4 


0.5 


20% 


1024 3 


K10-Q1 


55.2 


3.3 


0.83 


1.4 


1 


20% 


1024 3 



Karl et al. 


(2010 


I and 


Karl et al. 


(2011 



Masses refer to single galaxies in units of 10 Mq , unless stated otherwise. 



dark matter halo with concentration parameter c = 15, ex- 
ponential stellar and gaseous disks and a non-rotating Hern- 
quist bulge with masses [M ha i = 5.0 x 10 11 Mq, Af d isk,tot = 

4.1 x 10 10 Af ,M bu i gc = 1.4 x 10 10 M Q ], respectively. The 
disk is assumed to initially contain a fraction / g = 20% in 
gas (AW,* = 3.3 x 1O 1O M , M disk , gas = 8.3 x 1O 9 M ). 
Here we present a total of five runs, each with varying ef- 
ficiencies of the stellar feedback using a softened equation 
of state parameter with values of gEoS = (0, 0.01, 0.1, 0.5, 1), 
as detailed in Table [I] We use identical parameters for both 
progenitor disks except for the halo spin parameter A which 
directly influences the disk scale length of the gala xies (|Mo 
et al.|1998 1 and, hence, the length of the tidal tails (Springel 
fc White|1999| |. Our choice of A 40 38 = 0.10 and A4039 = 0.07 
results in disk scale lengths of rd,403s = 6.28 kpc and 
rd.4039 = 4. 12 kpc (see |Karl et al.|2010| for more details). 

At our standard resolution, each galaxy consists of a 
total of A/tot = 1.2 x 10 6 particles with 4 x 10 5 dark matter 
particles, 2 x 10 5 bulge particles, and a total of 6 x 10 5 
(stellar and gaseous) disk particles. The particle numbers 
are chosen such that we obtain a mass resolution of mb ary = 
6.9 x 10 4 Mq with a gravitational softening length ebary = 
35 pc in all baryonic components, while dark matter halo 
particles have a mass of 1.2 x 10 6 Mq with £dm = 150 pc. 
In the galactic disks, this choice yields 4.8 x 10 5 stellar and 

1.2 x 10 5 gas particles each. 

The initial interaction orbit is set identical to the one 
presented in |Karl et al.| ( |2010[ ). The progenitors are set 
on a mildly elliptic prograde orbit with initial separation 
r scp = r v ii = 168 kpc and pericentric distance r p = rd,4038 + 
fd,4039 = 10.4 kpc, where r v i r and r d are the virial radius 
and disk scale length, respectively. With disk inclinations 
*4038 = *4039 = 60° and pericentric arguments W4038 = 30° 
and W4039 = 60° this orbit has proven to result in a very good 
agreement with the observed large- and small-scale morphol- 
ogy and line-of-sight kinematics of the Antennae at the time 
of best match ( |Karl et al.|2010 |. 



Table 2. Model parameters of the radiative transfer calculations. 



2.2 Radiative transfer calculations 

In Section|3j we compare synthetic FIR maps and SEDs from 
our Antennae simulations with recent FIR observations by 
Klaas et al.| ( |2010 |, obtained with the Herschel- PACS instru- 
ment (see Section 2.3 1. We describe here the post-processing 
procedure applied to our simulations to construct the syn- 
thetic FIR maps and SEDs. 

The FIR emission in star forming galaxies is dominated 



Property 



RT parameter 



effective resolution 
constant gas-to-dust ratio 
dust model 
stellar emission model 
stellar IMF 
metallicity 
initial disk particles ages 
initial bulge particles ages 



1024 3 / 2048 3 
12 4 : 1 / 83 : 1 / 6 2 : 1 
|Draine| (|2003} 
BruzuaTSrUharlot[(|2003 1 



Halpeter 
7; 



(1955 



0.02 
- 6 Gyr 
5 - 9 Gyr 



by thermal emission from dust that has been heated by 
stellar radiation. Due to strong spatial fluctuations in both 
the dust density and the stellar emission, the dust tem- 
perature varies in a very complicated manner. To produce 
meaningful FIR maps from the simulations we model the 
three-dimensional radiative transfer (RT) separately with a 
dust RT code (|Juvela fc Padoan||2003"| |Juvela||2005| |Lunt-| 
tila fc Juvela|2012 1. Before starting the RT calculations, the 
SPH snapshot is gridded onto an adaptive mesh refinement 
(AMR) grid using the SPH smoothing kernel. The use of an 
adaptive grid enables good spatial resolution in the dense 
inner parts of the galaxy while keeping the total number 
of computational cells in the domain low. For most of the 
simulations described in this paper the size of the parent 
grid was chosen to be 80 kpc with an effective resolution 
of 1024 3 , resulting in a minimum cell size of ~ 78 pc. For 
simulations K10-Q0 and K10-Q0.1, we have performed ad- 
ditional tests with an effective grid resolution of 2048 3 and 
minimum cell size of ~ 39 pc without significantly changing 
the results. Details about the resolution of the RT calcu- 
lations are listed in Table Q] The total number of cells in 



the mesh refinement grid is typically ~ 5 x 10 . Interstellar 
dust is assumed to be uniformly mixed with the gas, and a 
single dust model is used throughout the galaxy. As a start- 
ing point , we use the st andard (Rv =3.1) Milky Way dust 
model of |Draine| (2003 Q This model implies a gas-to-dust 
mass ratio of 124:1. In addition, we test two models with a 
larger dust mass. In these latter two models, the properties 
of the dust grains are identical to the standard MW model, 
but the dust mass is multiplied by a factor of 1.5 and 2, 
yielding gas-to-dust mass ratios of ~ 83 : 1 and 62 : 1, re- 



http: / /www. astro. princeton.edu/>~draine/dust/dustmix. html 
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spectively. Stellar emission is modelled by assigning [Bruzual] 
& Chariot ( 2003 1 SEDs to all star particles according to their 



ages using a Salpeter stellar initial mass function (Salpeter 
1955). We assume Z — 0.02 ("solar abundance"; Anders 
& Grevesse 1989 1 for all stellar particles, i.e. disk, bulge, 



and newly formed stellar particles, which is in reasonable 
agreement with estimates from young star clusters in the 
Antennae (Bastian et al. 1:2006 2009). Stellar particles exist- 



ing at the start of the SPH simulation are assigned initial 
ages uniformly distributed between and 6 Gyr for the disk 
stars and between 5 and 9 Gyr for the bulge. Stellar parti- 
cles formed over the course of the simulation have their ages 
assigned corresponding to their formation time. The choice 
for the particular disk and bulge stellar age distributions has 
little influence on the FIR emission because the dust heat- 
ing is dominated by young stars, i.e., the stars formed in 
the simulations, depending on their specific star formation 
histories. 

After the dust distribution and radiation sources are set 
up, we calculate the dust equilibrium temperatures, assum- 
ing that the gas (and dust) distribution is smooth on the 
scale of individual grid cells, i.e. we do not use any sub-grid 
model for Hn regions (c.f. |Jonsson et aiT]2 010 Haywar d et al.| 
2011 ). All calculations are run iteratively to account for the 



dust self- absorption and heating until the dust temperatures 
have converged. Using the computed dust temperatures, we 
can construct SEDs in a wavelength range of 50/im - 850/im 
and spatially resolved maps in the 70 /im, 100 /jm, and 160 
/im passbands by integrating the radiative transfer equa- 
tion along the line-of-sight. The maps are then convolved 
to the resolution of the Herschel-PACS instrument at the 
corresponding wavelengths to simulate the observations by 



Klaas et al. (20101. The full- width-half-maximum (FWHM) 



of the point-spread function (PSF) is 5.5", 6.8", and 11.8" 
in the 70 pm, 100 /im, and, 160 /im bands, respectively. For 
comparison with the observations, we re-sized the simulated 
maps from the native RT resolution to the angular pixel sizes 
used in the observations (1.1", 1.4" and 2.1"), assuming a 
distance D — 30.8 Mpc to the Antennae. In addition, ran- 
dom Gaussian noise is added on top of the simulated data 
at a level corresponding to the noise levels in each spectral 
band, and we analyze the same spatial extent in the sim- 
ulated and observed maps. A summary of the RT model 
parameters is given in Table [2] 



source flux calibration overviewrl and the PACS photome- 
ter point spread function reporlrj). The derived fluxes are 
thus systematically lower by ~ 15% than those published 



Klaas et al. (20101. The final photometric uncertainties 



of the new maps are dominated by the current PACS abso- 
lute flux calibration accuracy, which is 3%, 3%, and 5% at 
70 fim, 100 fim, and 160 /im, respectively^. In the 160 /im 
filter, the [CII] 158 fim line contributes a few percent to 
the total flux as revealed by PACS spectra of the Anten- 
nae ( Gracia-Carpio et al. 20111. This means that the true 



160 fim continuum flux is likely lower by 1-2% than the one 
quoted in Table 4. 

Together with the 15 and 24 fim integral fluxes pub- 



lished in Klaas et al. (20101, 250 /jm, 350 fim, and 500 /mi 



photometry extracted from archival Herschel maps with the 
SPIRE instrument ( |Griffin et aT1|2010[ ), as well as 1.2 mm 
photometry from an IRAM-MAMBO map, the PACS 70 fim, 
100 fim, and 160 fim photometry is used to constrain the 
10 to 1000 fim spectral energy distribution of the Anten- 
nae galaxies as shown in Figure [l] This SED is decom- 
posed into three simple modified black bodies (see Figure 
[T|, which allows to apply proper colour correction factors 
for the PACS and SPIRE photometry. Furthermore, the sub- 
millimeter photometry enables an accurate prediction of the 
unobscured 850 fim flux which is used to estimate the total 



dust mass according to formula (6) in Klaas et al. (20011 



With a dust temperature of 23 K and an emissivity index of 
13—2, this gives a dust mass M dust = 7.1 x 10 7 M Q (for D = 
30.8 Mpc). 



Similar to the IRAS FIR flux definition (Helou et al 
1 1988} Lonsdale & Helo u] |1985[) and its e xtension to the 



170 fim ISOPHOT filter ( |Stickei et al.|20"00| ), we use a PACS 
FIR flux parameter for the wavelength range 59 to 206 fim 
(effective band width of the 3 PACS bands together) , which 
is valid for the temperature range 20 - 80 K and emissivity in- 
dices n = - 2, being the ranges of interest for galaxies. The 
conversion formula for the calculation of the PACS FIR flux 
parameter from the PACS 70, 100 and 160 fim monochro- 
matic flux densities will be presented in detail in |Klaas et] 
|al.| ( |2013"l in prep.). The total FIR luminosity is derived from 
this FIR flux parameter in analogy to the IRAS far-infrared 



luminosity calculation in Sanders & Mirabel (19961. This is 



the input for the calculation of observational IR star for- 



mation rates (SFR) according to formula (3) in Kennicutt 
fl998l ). 



2.3 Herschel-PACS observations of the Antennae 

As observational reference data we use far-infrared maps at 



70 /mi, 100 fim and 160 fim ( Klaas et al 



the PACS instrument (Poglitsch et al 



2010) obtained with 



2010) on-board the 



Herschel Space Observatory ( jPilbratt et al.||2010[). H ere we 



do not use the original maps published in Klaas et al. (2010), 



but re-reduced ones obtained in the context of a FIR study 
of a sample of nearby closely interacting galaxies ( |Klaas et| 



al.|2013 in prep.). The new processing is done with the pro- 
gramme package Scanamorphos ( Roussel|2012 I , resulting in 
an improved signal-to-noise (S/N) ratio and hence somewhat 
fainter levels in the flux limits (by factors of 1.84, 1.67 and 
1.97, respectively). Furthermore, more accurate calibration 
factors for extended emission are applied which were not 
yet available in early 2010 (see the PACS photometer point 



3 STAR FORMATION IN THE ANTENNAE 
SIMULATIONS 

3.1 Global star formation properties 

The top three panels of Figure [2] show colour-coded maps of 
the total gas surface densities, E gas , indicating differences in 
the central morphologies obtained for three Antennae sim- 
ulations providing a varying degree of pressure support in 
the star forming gas by adopting three different equation of 



2 http:/ /herschel. esac.esa.int/twiki/pub/Public/ 
PacsCalibrationWeb/pacs_bolo _fluxcal_report_vl.pdf 

3 http: / /herschel. esac.esa.int / twiki /pub/Public / 
PacsCalibrationWeb/bolopsL20.pdf 
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Table 3. Local star formation rates, gas masses, stellar mass formed within the last 15 Myr, and FIR fluxes in the galactic nuclei and the 
overlap region for the five Antennae simulations (K10-Q1, K10-Q0.5, K10-Q0.1, K10-Q0.01, and K10-Q0) with <j EoS = 1.0, g EoS = 0.5, 
ijEoS = 0.1, cjEoS = 0.01, and q Eo S = 0. 



Simulation 


Region 


SFRtMoyr" 1 ] 


M gas [xlO 9 M ] 


M*(< 15 Myr) [xlO 7 M ] 


/70//100//160 [JyQ 




NUC4038 


1.71 


1.60 


2.30 


4.59 / 6.54 / 4.12 


ijEoS = 1.0 


NUC4039 


4.08 


1.68 


6.30 


12.6 / 17.1 / 8.01 




Overlap 


0.49 


2.58 


1.03 


2.61 / 5.15 / 7.98 




NUC4038 


4.55 


1.49 


6.98 


10.3 / 14.3 / 6.63 


9EoS = 0.5 


Nuc 40 39 


2.85 


1.15 


4.72 


1U.4 / 12. 1 1 5. /4 




Overlap 


u. / D 


O.UO 


1 .00 






NUC403S 




U.oo 


1 A Q 


A AS2 / K OA I O SI 
4.40 / / Z.ol 


9EoS = 0.1 


NUC4039 


2.02 


0.10 


3.39 


9.91 / 11.8 / 5.65 




Overlap 


26.09 


3.67 


16.23 


22.5 / 39.4 / 35.3 




NUC4038 


0.49 


0.26 


0.75 


3.03 / 3.29 / 1.53 


q EoS = 0.01 


Nuc 40 39 


0.44 


0.14 


0.55 


3.67 / 3.76 / 1.59 




Overlap 


10.88 


3.09 


10.22 


29.7 / 39.3 / 27.6 




NUC4038 


0.99 


0.13 


1.48 


4.16 / 3.92 / 1.47 


ijEoS = 0.0 


NUC4039 


0.46 


0.15 


0.70 


3.31 / 3.44 / 1.35 




Overlap 


2.24 


2.52 


3.04 


8.44 / 11.4 / 10.1 




NUC4038 


1.57 ±0.04 


3.65 




3.91 ± 0.12 / 6.73 ± 0.20 / 8.94 ± 0.45 


Observation^ 


Nuc 40 39 


0.72 ± 0.02 


1.34 




1.91 ± 0.06 / 3.53 ± 0.11 / 3.40 ± 0.17 




Overlap 


6.93 ±0.27 


6.16 




21.6 ±0.4 / 31.4 ±0.7 / 30.1 ± 1.1 



a Fluxes at 70/^m, lOO^tm, and 160/^m, given in Jy, for models with gas-to-dust rati o 62 : 1 and the ob servations (see Section 2.3 1 
^ SFRs are calculated from the total infrared luminosity, L( 50 _ iooo)Mm; according to 



Kennicutt 



1 1998 1, equation (3). FIR flux densities, 



gas masses, and SFRs for the overlap region are estimated as a lower limit by summing up the emission from knots Kl — K4 (for D 
30.8 Mpc). See|Klaas et al.|j2010| for a definition of the different FIR knots. Observed gas masses denote molecular gas masses from 
Klaas et al.|2010 



Wilson et al. 



(20001, see 



Table 2. 



state parameters (see Sectionj2|, q^oS = 1-0, q Eo s = 0.1, and 
gEoS = 0.01 (from left to right). Each gas surface density 
map is overlaid with contours of the SFR surface density 
Ssfr- For all three parameters, we find a similar large-scale 
morphology: the nuclei of the progenitor disks are still well 
separated, but connected by a ridge of high density gas. In 
the case of very efficient feedback (?EoS = 1-0; K10-Q1.0), a 
high-density overlap region forms, but the highest densities 
are found in the galactic nuclei. In the simulations K10-Q0.1 
and K10-Q0.01 with less vigorous feedback (gEoS — 0.1 
and qeoS = 0.01), the high surface density regions tend 
to become tighter and more pronounced in the overlap 
(middle and right top panels). As a consequence we also 
find an increasing number of low-density ("void") regions 
inside the disks for lower values of qeoS ^ 0.1 due to 
the decreasing pressure support in the star forming gas. 
The star formation contours directly support this picture. 
Lowering the effective pressure, and hence the stellar 
feedback, in star forming regions leads to more confined, 
localized regions of intense star formation activity. In the 
bottom panels of Figure [2] we show the corresponding 
contours of the gas surface density and overplot all stellar 
particles (with individual masses of m s tar = 6.9 x 10 4 M©), 
colour-coded according to their ages, having formed in the 
last r < 15 Myr (blue), 15 Myr < r < 50 Myr (green), and 



50 Myr < r < 100 Myr (red J] This period corresponds to 
the time span of the interaction-induced rise in the star 
formation activity during the recent second encounter. 
The youngest stars (blue) form predominantly in regions 
of currently high gas densities, i.e. in the centres, the 
overlap region, and in the arc-like features along the disks, 
very similar to observations of the youngest clusters in 



the Antennae (Whitmore et al. 1999 2010 Wang et al 



2004b 1. This qualitatively confirms our earlier results in 
Karl et aL| ( |2010| , indicating that our conclusions are robust 



to altering the parameters of the adopted star formation 
feedback algorithm. However, Figure [2] also shows that 
the star formation in the galaxy centres is much more 
pronounced in run K10-Q1.0 with qeoS = 1.0 (bottom 
left panel), while the overlap starburst gains more relative 
importance for runs K10-Q0.1 (qeoS = 0.1) and K10-Q0.01 
(<?EoS = 0.01; see also below). Furthermore, young stars 
(r < 50 Myr) tend to be more concentrated in the overlap 
region and the well-defined arcs in the remnant spirals, 
while older stars (50 Myr < r < 100 Myr) are dispersed 
more evenly throughout the disks. This captures a similar 
trend observed in the Antennae (e.g. |Zhang et al.|[200l| 
their Figure 2). Therefore we conclude that, within a range 
of feedback efficiencies (?EoS < 0.1, the overlap starburst 



4 The corresponding stellar masses formed in the last r < 15 Myr 
are given in Table [3] for different regions in the merger. 
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Table 4. Global star formation rates, gas masses, and FIR fluxes, 
measured over the entire system, for all Antennae simulations 
(K10-Q1, K10-Q0.5, K10-Q0.1, K10-Q0.01, and K10-Q0) with 
varying feedback parameters: QeoS = 1-0, 5EoS = 0.5, QeoS = 0-1) 
9EoS = 0.01, and q E oS = 0.0. 



Simulation 


SFF0 




/tot, 70 / /tot, 10C)//tot, 16C 


b 




QEoS = 1-0 


7.16 


13.2 


31.0 / 49.6 / 46.8 


QEoS = 0-5 


8.74 


12.2 


34.0 / 51.7 / 45.1 


<3EoS = 0.1 


30.1 


11.8 


48.5 / 76.8 / 66.4 


<? EoS = 0.01 


14.3 


9.28 


49.2 / 65.1 / 46.7 


9EoS = 0.0 


5.36 


8.38 


26.8 / 34.0 / 26.7 


Observation^] 


19.1 ± 1.4 


~ 9.6 


56.5 / 89.4 / 86.1 



a Units for the numbers given in columns 2-4 are in [Mq yr 1 ], 
[xlO 9 M0], and [Jy], respectively. 

^ Total fluxes at 70^tm, 100/xm, and 160/^m, given in Jy, for models 
with gas-to-dust ratio 62 : 1. 

c For the observed SFRs and FIR fluxes we give values re- 
analyzed for this paper (see Section |2.3| Klaas ct al. 2 013| in 
prep.). Total gas masses in the galactic discs for H2 and Hi are 
adopted from |Gao et al.| ( |2001| ) and |Hibbard et"aT] ( |2001| > , assum- 
ing a distance of D = 30.8 Mpc and a CO-to-H2 conversion fac tor 
a co = 0.8 Mq pc~ 2 (K kms -1 )" 1 1 Downes fc Solomon j 1998 \ . 



seems to be driven by the most recent episode of star 
formation in the Antennae, which was induced by the 
second close encounter, while for values approaching the 
full stellar feedback (§EoS > 0.5) the overlap starburst is 
significantly suppressed. 

Thanks to the very good spatial resolution of the 
Herschel-PACS instrument in the FIR, the previously de- 
termined IR SEDs of individual star forming knots (see e.g. 
|Brandl eF al. 2009) could be augmented by additional flux 
estimates at 70pm, 100 /im, and 160 /im ( Klaas et al.|2010 



Klaas et al.|2013 in prep.). The improved analysis described 
in Section |2.3| allows for a precise determination of the to- 
tal FIR luminosities and SFRs in individual knots yielding 
integrated SFRs of SFR 4038 = (1.57 ± 0.04) Mq yr -1 and 
gFR 4039 = (0 72 ± o.02) Mq yr" 1 in the nuclei of NGC 4038 
and NGC 4039, and SFR° vL > 6.93 M yr -1 , when taking 
the sum over all knots in the overlap region as a lower limit 
(i.e. knots "Kl" to "K4"; see Figure 3 in |Klaas et "aI][20To| 
for a definition). We compare these estimates to our sim- 
ulated SFRs in the galactic nuclei of NGC 4038 and NGC 
4039 (defined as the central region of each nucleus with a 
radius of lkpc), and the overlap region (chosen according 
to the central morphology of the simulations). Details are 
provided in Table |3j where we compare the total SFRs, the 
gas mass, the stellar mass formed in the last 15 Myr, and 
the FIR fluxes in the three PACS bands at 70^m, 100/im, 
and 160/im for the different specific regions in each sim- 
ulation. In the case of efficient feedback (<?EoS = 1.0 and 
9EoS = 0.5), the simulations shows stronger star forma- 
tion activity in the two nuclei compared to the overlap re- 
gion. The simulations with low stellar feedback (§eoS = 0.1 



and gEoS = 0.01), develop a very vigorous overlap starburst 
with a SFR of ~ 26M yr -1 and ~ llM yr -1 , which is 
higher than the observed lower limit of > 6.93 Mq yr -1 (by 
a factor of ~ 1.5 — 4) together with a very low ratio of 
SFRnudci/SFRovoriap = 0.11 and 0.09 compared to < 0.33 in 
the Herschel-PACS observations. Simulation K10-Q0 with- 
out stellar feedback has the lowest gas masses and SFRs of 
all simulations due to a very efficient star formation and gas 
consumption earlier in the simulation (see below). The to- 
tal instantaneous SFRs in the five simulations show a large 
scatter within 5.4M yr -1 < SFR < 3OM yr -1 , as mea- 
sured from the simulated SPH particles at the time of best 
match (see Table B. These values are in broad quantita- 
tive agreement with t he range of observ ed values between 
~ 5-20 M yr -1 (e.g. |Zhang et al.|200l"| and with the value 
of SFR = (19.1 ± 1.4)M yr -1 we derive from the total IR 
luminosity in this paper (see Table Q. Therefore, without 
attempting a fine-tuned match, we conclude that adopting 
a rather inefficient parametrization for the thermal stellar 
feedback, such as in simulations K10-Q0.1 and K10-Q0.01, 
tends to give the most consistent results with the observed 
spatial distribution of star forming regions in the Antennae 
system. 



3.2 Reproducing the overlap starburst 

In a first attempt to explain the large scatter in the star for- 
mation properties, we show the evolution of the gas masses 
for three simulations employing different parameters for the 
stellar feedback model in Figure [3] The gas is separated into 
non-star-forming (blue) and star-forming gas phases (or- 
ange) based on the density criterion n > n cr i t = 0.128cm -3 . 
All three simulations show a steady decrease in their gas 
mass (upper solid line) due to consumption by star forma- 
tion from the very start of the simulations, corresponding 
to a gradual build-up of the new stellar component (blue 
solid line). However, due to the lower pressure support in 
the effective multi-phase model in the K10-Q0.01 simula- 
tion (bottom panel), a fraction of 12.5% of the initial gas 
mass has already been turned into stars by the time of 
the first pericentre (blue dots), whereas the K10-Q1.0 sim- 
ulation, with a more efficient stellar feedback (top panel), 
has only consumed about 7.2% at the same time. The K10- 
Q0.1 case is intermediate with 10.4%. The total fraction of 
non-star-forming gas (blue) stays constant to within ~ 10% 
around a level of 0.3 during the whole course of the simu- 
lated time span (t < 2.0 Gyr), save some little 'dips' during 
episodes of induced star-formation associated with the first 
and the second pericentre passages indicated by the blue 
(t = 0.66 - 0.68 Gyr) and green (t « 1.20 Gyr) dots in Fig- 
ure H The more dramatic evolution can be seen in the star 
forming gas (orange, n > n cr it = 0.128 cm -3 ) which gets ef- 
ficiently depleted after the first pericentre in the K10-Q0.01 
simulation with weak feedback (bottom panel), having con- 
sumed more than 41% of the initial gas mass in the form 
of star forming gas at the time of second pericentre. In con- 
trast, in the full-feedback run (K10-Q1.0) there is still about 
81% of the initial gas mass remaining after the second en- 
counter, which leads to a very rapid consumption of 40% 
of the gas within the first ~ 250 Myr after the final merger 
(red dots). The key difference is that the low feedback sim- 
ulation consumes the nuclear gas reservoir before the best 
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Time [Gyr] 

Figure 3. Time evolution of the mass in gas (upper black solid line) and newly formed stars (blue line) normalized to the initial gas mass 
in three Antennae simulations with different parameters for the softened equation of state. Top panel: full feedback simulation K10-Q1 
with qEoS = 1-0; middle panel: intermediate feedback simulation K10-Q0.1 with qeoS = 0.1; and bottom panel: the K10-Q0.01 simulation 
with very weak feedback, <?EoS = 0.01. In addition, we distinguish between the dense star forming gas (n > n cl .; t = 0.128 cm -3 ) treated 
by the effective multi-phase model (orange) and non-star forming gas (blue). The blue, green and red dots indicate the times of first 
pericentre, second pericentre, and the time of best match in each simulation. 
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Figure 4. Star formation histories for the five Antennae simulations with varying feedback parameters: QeoS = 1-0 (green), <?EoS = 0-5 
(purple), QEoS = 0.1 (blue), 5eoS = 0.01 (red), and 5EoS = 0.0 (orange). The arrow gives the time of first close encounter, while the bar 
indicates the time between second close encounter and final merger. The time of best match is denoted by the coloured dots for each 
simulation. 



match in contrast with the full feedback simulation. In Ta- 
ble [3] we find that, with decreasing gEoS, the nuclear gas 
masses go down drastically — from nuclear gas masses of 
M gas = 1.6 x 10 9 M Q for q Eo s = 1.0 to a factor of ~ 10 
lower in the other two runs. The gas masses in the overlap 
regions, on the other hand, show a less pronounced depen- 
dence with the chosen effective stellar feedback parameter, 
staying constant to within ~ 40%, while the SFRs increase 
significantly for the two lower-feedback simulations as dis- 
cussed above. 

A complementary view of this picture is offered by the 
star formation histories in the Antennae simulations (Fig- 
ure Q. The isothermal and quasi-isothermal simulations 
(K10-Q0 and K10-Q0.01), show quite similar initial SFRs 
of ~ 3.6M yr _1 and ~ 3.2M yr _1 . The SFRs for the 
simulations with efficient feedback and qeoS 0.5 (K10- 
Q0.5 and K10-Q1.0) are smaller by a factor of about ~ 1.5 
from the very start, while run K10-Q0.1 is intermediate 
with SFR(t = 0) ~ 2.6M©yr~ 1 . After the first pericentre 
(black arrow) runs K10-Q0 and K10-Q0.01 both experience 
a significant rise in their SFRs by more than a factor of 
five, while the other three simulations show a more grad- 
ual increase. The low-feedback simulations show high total 
SFRs of ~ 30 M Q yr _1 (K10-Q0.1) and ~ 14 Mq yr" 1 (K10- 
Q0.01) at the time of best match (see Table |4J, powered 
by the short-lived overlap starburst shortly after the second 
pericentre. The full-feedback (K10-Q1.0) and intermediate- 
feedback (K10-Q0.5) simulations, on the other hand, build 
up very powerful nuclear starbursts after the final merger 
with a maximum total SFR of SFR max ~ 14OM yr _1 
and SFRs at the best match of only ~ 7M0yr _1 and 
~ 9M yr _1 (Table Ek, respectively. An equally power- 
ful starburst is missing in the case of weaker stellar feed- 



back with SFR max = 68.4 Mqyv- 1 for the K10-Q0.1, and 
SFR max = 4O.8M yr _1 for the K10-Q0.01 simulations. 
Having consumed about half of its gas before the second 
pericentre, simulation K10-Q0 has the lowest SFR of all 
simulations during the time between second pericentre and 
final merger with a SFR ~ 5.4Af yr at the best match. 
Leaving aside simulation K10-Q0 without stellar feedback, 
we find a clear anti-correlation between the height of the 
maximum SFRs after the first and second pericentres in the 
sense that simulations having higher star formation peaks af- 
ter the first pericentre have correspondingly lower star for- 
mation peaks after final coalescence and vice versa. Note 
also that this trend, interestingly, yields very similar gas 
fractions in these simulations after final coalescence of the 
galaxies (see Figure [31, differing by less than ~ 5% at 2 
Gyr. It is, however, not directly reflected in the SFRs at the 
time of best match, which we find shortly (~ 30 — 50 Myr) 
after the second encounter but before the time of maximum 
star formation activity (see Figure Q. This is firstly due to 
a more rapid gas consumption during the first pericentre in 
the weaker feedback runs and secondly due to a weaker pres- 
sure support in the star forming gas phase in the K10-Q0.1 
and K10-Q0.01 simulations making the gas in the overlap re- 
gion more susceptible to local shock-induced fragmentation 
and subsequent star formation, as discussed above. 

Overall, the star formation histories are complex and 
there is a large scatter in the SFRs at the time of best match 
(being defined as the time when simulated and observed 
morphology and line-of-sight kinematics are most similar), 
depending on the parameter choice of (/Eos for the stellar 
feedback model. This, e.g., leads to very similar SFRs at the 
time of best match in the full-feedback simulation K10-Q1, 
employing a vigorous feedback, and the isothermal Simula- 
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tion K10-Q0 with no feedback at all. Therefore, global prop- 
erties like the total SFRs and the total gas masses alone seem 
to be a degenerate diagnostic when comparing to observa- 
tions. 



3.3 Comparison with observed CO maps 

As a first reliability check of the hydrodynamical simula- 
tions, we construct a simple CO- map derived from the simu- 
lated star formation rates for comparison of the simulations 
to this widely used molecular gas tracer. Hereby, we first 
invoke the amount of molecular H2 from the SFRs of the 
simulated gas particles by using an "inverse" |Bigiel et al.| 
( 2008 ) relation of the form 



J H 2 



10- 



A 



l/N 



M pc" 



(1) 



where £jj 2 and Esfr are the surface densities in 
H2 and the SFR, respectively. The latter is given in 
units of M© yr _1 kpc -2 , A is the SFR surface density 
at a fiducial gas surface density of 10 Mq pc -2 with 
Iog(A[M yr" 1 kpc" 2 ]) = -2.06, and N = 0.96 is a power- 
law index near unity (see B igiel et al.|2008 Equation 2 and 
Table 2). To calculate the simulated SFR surface densities, 
we bin the gas particles on a two-dimensional grid with 
pixel sizes corresponding to 1" x 1" at the assumed dis- 
tance of D = 30.8 Mpc. Then, we relate the H2 masses (per 
pixel) back to integrated CO intensities (per pixel) by adopt- 
ing the conversion M mo i = 1.61 x 10 4 (D/M pc ) 2 Sco M Q 
(Wilson & Scovillell 19901 1 Wilson et al.|2000b, where Sco is 



given in Jy tons . Finally, we convert the map to units of 
Jy kms _1 beam _1 by multiplying the whole map with the 
beam area in pixels. The effective beam area is given as 
1.133 x (3". 15 x 4" .91) (FWHM) (Wilson et al.|2000| > where 
the factor 1.133 is a correction for the Gaussian shape of 
the bearr{^] We show the synthetic CO map for simulation 
K10-Q0.01 in the left panel of Figure [5] together with the 
observed integrated CO map to the right (taken from Fig- 
ure 1 in Wilson et al.|2000 |. Here, we have applied a Gaus- 
sian filter to the simulated CO map, corresponding to the 
FWHM sizes of the synthesized beam (see above). We use 
the same contour range as the observations. However, since 
the simulated CO emission is more concentrated (see also 
Section [4] for a discussion) we put an additional contour at 
0.1 Jy beam -1 kms -1 to account for the lack of extended 
emission in the simulation. The three most luminous peaks 
in the simulation actually show higher peak values (by up 
to a factor of five) than the highest contour level in the ob- 
servation. These can be seen as the blank spots within the 
three peaks in the overlap region. In the synthetic CO map 
(left panel) we find the CO emission to be mostly confined 
to the central regions of the merger and a number of dis- 
tinct high-density peaks in the southern overlap and near 
the two nuclei. To a lesser extent, i.e on the additionally 



5 One should bear in mind, however, that these maps still fun- 
damentally represent the simulated SFRs and no CO(1-0) emis- 
sion, which we do not model in the simulations presented here, 
even if we will call the converted simulated SFR maps 'CO-maps' 
(instead of 'SFR-to-CO-maps') for simplicity for the rest of this 
paper. 



added, lowest level contour, we also find emission from two 
arcs along the southern and northern remnant arms, plus 
one prominent peak in the southern arm. The observations 
exhibit a generally similar picture (right panel in Figure j5J). 
The strongest emission peaks from the CO (1-0) transition 
are also detected in the two galactic nuclei and the overlap 
region. However, as noted above, there is some discrepancy 
in the more extended emission, in the sense that the north- 
western spiral arc in NGC 4038 shows about 10 times higher 
emission than the simulation, and, there is no emission de- 
tected in the south-western region of the southern disk as 
found at a low level and in one distinct peak in the simu- 
lation. The apparent CO excess in the southern simulated 
disk is qualitatively consistent with the high total gas sur- 
face densities in Figure [2] as discussed above and noted in 
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4 COMPARISON WITH FAR-INFRARED 
OBSERVATIONS 

In Figure [6] we compare Herschel-PACS data (top row; 
see Section |2.3[ ) in three photometric wavebands centred 
on 70/im, 100/xm and 160/im with the simulated FIR RT 
maps of the K10-Q0.5, K10-Q0.1, and K10-Q0.01 simula- 



tions (bottom three rows; see Section 2.2 1. Being an a priori 



un-constrained parameter in the RT model (see Section 2.2 1, 
we tested for variations in the gas-to-dust ratio (not shown) 
but found only negligible changes in the FIR morphology 
for gas-to-dust ratios between 62 : 1 and 124 : 1 in all three 
FIR bands. Therefore, we only show the simulated maps 
adopting a gas-to-dust ratio of 62 : 1, which generally gives 
the best match to the observed SED in the Antennae (see 
below). We choose a lower cut in the surface brightness de- 
termined by the mean noise level corrected for correlated 
noise above the mean background in each band. The upper 
limits in the colour scheme are set by the maximum surface 
brightness levels of the observations. 

Emission in the FIR bands reveals the sites of recent 
star formation that are still embedded in the dusty molecu- 
lar clouds from which they have formed. The simulated FIR 
properties vary significantly depending on the adopted ef- 
ficiency of the sub-grid stellar feedback model. Simulations 
K10-Q0.5 and K10-Q0.1 (two bottom rows) show only a 
small number of emission peaks, most dominant either in the 
two nuclei (K10-Q0.5) or in the overlap region (K10-Q0.1) 
in all three simulated bands. Simulation K10-Q0.01 (second 
row of Figure [6] on the other hand, shows a well-defined 
pattern of sequential, confined star forming knots along the 
overlap region very similar to the observations (top row). 
However, there are also two additional star-forming knots 
near the northern nucleus which are not detected at a sim- 
ilar level in the observations. As seen in Tables [3] and [4] 
more than > 50% of the FIR emission in each band origi- 
nates from the overlap region in simulations K10-Q0.01 and 
K10-Q0.1, compared to only 6 - 14% (12 - 30%) from the 
two nuclei in NGC 4038 and NGC 4039 combined in K10- 
Q0.01 (K10-Q0.1). A similar percentage (10-15 %) is ob- 
served in the combined fluxes of the two nuclei in all three 
Herschel-PACS bands. The combined emission from the ob- 
served knots in the overlap regions gives a lower bound for 
the total emission from the overlap region, yielding a total of 



Antennae simulations compared with Herschel-P^ CS 11 




_5 5 a(J2000) 

x 1 [kpc] 



Figure 5. Left panel: Synthetic CO map of the central regions of the Antennae from the K10-Q0.01 simulation. The integrated CO 
masses are obtained in a two-step process. First, the H2 surface densities are obtained from the SFR surface densities of the SPH particles 
using an "inverse" |Bigiel et al. | |2008| l relation. Then, adopting a conversion motivated by the observations l |Wilson et al,|2000| |, we map 
H2 masses back to CO integrated intensities (see text for details). Right panel: CO(l-O) integrated intensity map, adapted from Figure 
1 in |Wilson et al.] | |2000[ | . The two galactic nuclei and five super-giant molecular clouds in the overlap region are indicated. The filled 
black oval in the top left corner gives the size of the synthesized beam (3". 15 x 4". 91) in the observations. 



~ 35 — 39% for the three bands. We find a very good agree- 
ment between the observed FIR surface brightness maps and 
those of simulation K10-Q0.01 

Apart from differences in the spatial distribution of the 
simulated and observed FIR emission due to the specific 
adopted stellar feedback parameters, there are also some 
characteristic differences with respect to the observed FIR 
maps which are generic to all the simulations. These are 
mostly caused by slight differences in the underlying overall 
projected gas morphologies in the simulations, as discussed 
in Section [3] (see Figure |2l. Most notably, the simulated RT 
maps in the bottom three rows in Figure[6]mdicate some low- 
level emission from the southern disk (NGC 4039) which is 
not observed (top row). This makes the observed disks look 
slightly more tilted to the north-west than the simulated 
ones. Furthermore the dominant emission peak is found in 
the northern part of the overlap region in the simulations 
which is slightly at odds with the highest emission being 
observed in knots Kl and K2, at the southern end of the 
overlap region in the Antennae (Figure 3 
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20101. Comparing the global distribution of FIR emission, 



we do not find equally strong spatially extended emission in 
the simulated progenitor disks. This may be best seen in the 
northern spiral, e.g. comparing the arc-like feature which is 
clearly detected in the 100 (im and 160 /im Herschel maps 
(top middle and left panel) to the very localized emission 
in the simulations in these bands (bottom three middle and 
left panels). This may partly be caused by different large- 
scale gas morphologies, as discussed above, and, to some ex- 
tent, by resolution effects and/or restrictions in the adopted 
sub-grid model for the stellar feedback, leading to a spa- 
tially confined star formation in high-density regions, which 
in turn is seen in the simulated FIR maps. 

To highlight the striking differences in the simulated 
FIR properties when adopting different feedback parameters 



we present in Figure [7] the two most extreme cases from our 
sample, the full-feedback simulation K10-Q1 with ^eoS = 1 
and the isothermal simulation K10-Q0 without any stellar 
feedback, i.e. gEoS = 0, both with a gas-to-dust ratio of 
62 : 1. In the K10-Q1 simulation, emission from the overlap 
region is feeble, especially at the shorter wavelengths, which 
reveals a prevalence of cold dust in this region. The ob- 
served overlap starburst is missing entirely in the synthetic 
FIR maps. In simulation K10-Q0 a large number of spatially 
confined peaks of emission develops in the three FIR bands 
and the emission shows a secondary peak at the southern 
end of the overlap region. Due to the high prior gas con- 
sumption and a resulting low SFR of ~ 5.4M Q yr _1 at the 
time of best match, however, the surface brightness is much 
lower compared to the best-matching simulation K10-Q0.01 
(as discussed in Section I3J. 

Figure [8] shows the simulated SEDs of the three An- 
tennae simulations with qeoS = [1.0, 0.1, 0.01] in comparison 
to observations of the FIR fluxes in the Antennae at 70/im, 
100/im, and 160/im, given as the filled dots with error bars, 
and a model of three combined modified black body spectra 
fitted to IR data points within 10— 1000/im, given as the thin 
black line (see Figure [I] and Section 2.3 1. The different lines 
for each simulation probe varying gas-to-dust models with a 
standard Milky Way value of 124 : 1 (dotted lines) and lower 
values employing successively more dust with ratios of 83 : 1 
(dashed lines) and 62 : 1 (solid lines). Increasing the total 
amount of dust in the RT simulations leads to increasingly 
lower dust temperatures seen as shifts in the peak emission 
towards longer wavelengths. The observed spectral shape is 
best reproduced with a low gas-to-dust ratio of 62 : 1 (solid 
lines). In particular, this choice leads to a steep slope be- 
tween 70/xm and 100/xm, similar to the observed data points, 
with a ratio /loopm/iVo^m = 89.4/56.5 « 1.6, as well as 
the flat and subsequently declining part near the peak of 
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Figure 6. Synthetic 3D radiative transfer maps of the central regions of the Antennae at 70 /im (left column), 100 /im (middle column), 
and 160 /im (right column) for simulations K10-Q0.01, K10-Q0.1 and K10-Q0.5 with a gas-to-dust ratio 62 : 1 (bottom three rows). The 
simulate d maps are compared to Herschel-PACS observations in the same bands (top row; see Section [273] |KIaas et al.|2010| [KTaas et| 
|al.|2013| in prep.). Surface brightness is given in logarithmic units of Jy/arcsec 2 . 
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Figure 7. Synthetic 3D radiative transfer maps of the central regions of the Antennae at 70 /im (left column), 100 /im (middle column), 
and 160 /im (right column) with a gas-to-dust ratio 62 : 1, highlighting the large scatter in FIR properties when adopting the two extreme 
parameter choices for the stellar feedback model in our sample: the simulation without stellar feedback K10-Q0 with (JeoS = 0-0 (t°P) 
and the full-feedback simulation K10-Q1.0 with qeoS = 1-0 (bottom). Surface brightnesses are given in logarithmic units of Jy/arcsec 2 . 



the spectrum between 100pm and 160pm. Only for wave- 
lengths below < 60pm the simulated SEDs seem to drop 
a bit too fast compared to the fitted modified black body 
spectra. This is in part due to the fact that, at these and 
shorter wavelengths we start to miss some significant contri- 
butions to the total flux density from stochastically heated 
small dust grains, or, in the mid-infrared, from polycyclic 
aromatic hydrocarbons (PAHs). 

More important for the SEDs of the simulated data, 
however, is the choice of the adopted stellar feedback pa- 
rameter qeoS- While the adopted gas-to-dust ratio changes 
mostly the spectral shape, the parameter for the stellar feed- 
back changes significantly both the spectral shape and the 
total flux of the simulated SEDs due to differences in geome- 
try, total gas mass, and SFRs between the simulations. Sim- 
ulation K10-Q0.01 has the smallest amount of gas (and dust) 
in the actively star forming regions (see Table [3| combined 
with a relatively high SFR at the time of best match, yield- 
ing a strong emission from hot dust at shorter wavelengths. 
Simulations K10-Q0.1 and K10-Q1.0, however, show a spec- 
tral shape peaking at longer wavelengths due to their larger 
reservoirs of gas (and dustj^j The best agreement to the ob- 



served data points both in terms of the total flux level as 
well as the spectral shape is reached for simulation K10-Q0.1 
(i.e. with a low feedback efficiency of gEoS = 0.1) and a gas- 
to-dust ratio of 62 : 1 (red solid line in Figure [8]). The SEDs 
of simulations K10-Q0.01 and K10-Q0.1 peak at a value of 
~ 65 Jy at 99pm and ~ 81 Jy at 116pm (compare also with 
Figure [(Sj|, respectively, similar to the peak value of 97.4 Jy 
at 118pm in the modified black body spectrum (thin black 
line) . 



5 SUMMARY & DISCUSSION 

In this paper, we analyze the global and spatially resolved 
star formation properties in a suite of hydro-dynamical sim- 
ulations of the Antennae galaxies in order to understand the 
origin of the observed central overlap starburst. Our aim is 
to assess the impact of variations in the applied thermal stel- 
lar feedback on observable simulated star formation proper- 
ties at the time of best match, in the spirit of a study on the 



Mice galaxies by Barnes (20041 who investigated variations 



6 Note here that average dust temperatures depend much more 
strongly on the total dust mass than on the SFR, or the bolo- 



metric luminosity L\, a \, causing the SED to appear "warmer" 
smaller dust masses (e.g. |Hayward et al.|2011[ ). 



for 
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Figure 8. Comparison of the total spectral energy distribution of radiative transfer calculations for the three different Antennae 
simulations K10-Q1, K10-Q0.1, and K10-Q0.01 with observations. Dotted lines correspond to SEDs calculated assuming a standard 
Milky Way gas-to-dust ratio of 124:1, while for the dashed and solid lines the ratio is reduced to a factor of 83 : 1 and 62 : 1, respectively. 
The small circles give the total fluxes obtained with Herschel-PACS observations at 70 /im, 100 (im, and 160 /im (see Figure[l]and Section 
|i01|Klaas et al.|2010||Klaas et al.|2013[ in prep.). The total (intrinsic+systcmatic) photometric error, is given as bars, where we estimated 
the intrinsic photometric error as ~ [50, 50, 90] mjy and the systemic error as ~ [3%, 3%, 5%] at 70 /im, 100 fim, and 160 /im, respectively. 
The thin black line gives a theoretical fit to the observed IR, SED using a three-component modified black body model (see Section|2.3|. 



due to the applied star formation law. Using a state-of-the- 
art dust radiative transfer code ( Lunttila fc Juvela|2012 



Section 2.2 1 we directly compare synthetic FIR maps and 



SEDs with re-analyzed Herschel-PACS observations ( Klaas 
et al.|20T0| see Section [23] Klaas et al.|2013 in prep.). This 



is the first time that the FIR properties of a specific local 
merger system are compared in quantitative detail to the 
simulated FIR properties of a dynamical merger model. 

The distributed extra-nuclear starburst in the Antennae 
galaxies is plausibly a natural consequence of induced star 
formation in high density regions in the overlapping galac- 
tic disks after the recent second close encounter. This is the 
only phase in the simulations when a gas-rich, star-bursting 
overlap region forms and when we find spatially confined ar- 
eas of high star formation activity and of similar extent as 



observed in the overlap region in the Antennae (e.g. [Mirabel 
et al.|1998||Zhang et al.|2001||WiIson et al.|2000||Wang et ~ 
2004b||Zhang et al.|2010[ ). We confirm that this result holds 



for a range of adopted stellar feedback efficiencies of the sub- 
grid star formation model (see e.g. Figure Ejb. However, we 
find a trend that the overlap star formation activity is sig- 
nificantly suppressed for values approaching the full stellar 
feedback simulation ((/BoS < 1; Springel fc Hernquist|2003 1. 
Best qualitative agreement with observations is obtained us- 
ing a rather inefficient parametrization, 0.01 < §eoS < 0.1, 
for the thermal stellar feedback. This is in agreement with 



recent results by |Hopkins fc Quataert (20101 who advocated 
a range of qeoS = [0.125,0.3] by comparing model predic- 
tions and simulations for different values of qEoS to the gas 
properties in observed star forming systems. 

The star formation histories in the different simulations 
are complex and there is a large scatter in the global ob- 
servable properties, e.g. the total gas mass, star formation 
rate and FIR fluxes, at the time of best match (see Table [4] 
Figures [3] and [4| , even though they only differ in the cho- 
sen parameter for the thermal stellar feedback. We find that 
additional information provided by a spatially resolved anal- 
ysis of these properties is needed to break degeneracies in 
the model assumptions. 

Synthetic FIR maps at 70/xm, 100/im, and 160/im (Fig- 
ures [6]&[7f show very different FIR surface brightness prop- 
erties for varying efficiencies of the thermal stellar feedback. 
We find that, for a rather inefficient stellar feedback for- 
malism with 0.01 < qeoS < 0.1, the simulated maps are 
in good agreement with the spatially resolved FIR obser- 
vations. This quantitatively supports results from a simple 
conversion of the simulated local star formation to molecular 
hydrogen surface densities that compares well with the ob- 
served integrated CO intensity map (Figure [5|. We propose 
that the reason for this is that the simulations employing 
low feedback efficiencies encompass a sweet spot between a 
moderate previous gas consumption history together with 
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a higher susceptibility to induced star formation after the 
second close encounter caused by the weak thermal support 
against local gas collapse. 

Probing synthetic SEDs in the FIR using different gas- 
to-dust ratios, we find that a low gas-to-dust ratio of 62:1 
yields spectral shapes of the simulated SEDs that are in 
good agreement with a simple modified black body spec- 
trum fitted to the observed data points; however, with sig- 
nificantly lower flux densities in the simulated SEDs (Fig- 
ure)^. For simulations K10-Q1.0, K10-Q0.1, and K10-Q0.01 
with gas-to-dust ratios of 62 : 1, the average flux decrements 
are 43.9%, 18.6%, and 31.7% within 50 - 280/^m. The syn- 
thetic far-infrared SED of the best-fitting simulation, K10- 
Q0.1 (i.e. with gEoS = 0.1), peaks at a value of ~ 81Jy 
at 116/im, similar to the peak value of 97.4 Jy of the fitted 
modified black body spectrum at 118/xm. A comparison with 
the observed L 25 o-SFR relation (Hardcastle et al.||2013| c.f. 



formulae on page 13; their Figure 9) reveals that the SFRs 
estimated from our modelled flux densities at 250/im seem 
to be systematically lower (by an average of 0.4 dex) than 
the SFRs in the hydrodynamic simulations. About half of 
the simulations lie within the observed uncertainty in the 
SFRs, which is given as ±0.5 dex ( jHardcastle et aL||2013[ 
their Figure 9). Again we find that the agreement is best for 
simulations with low gas-to-dust ratios of 62 : 1. 

Is the seemingly favoured low gas-to-dust ratio of 62:1 
realistic? Assuming all the metals to be locked up in the dust 
grains and none left in the gas phase, this value corresponds 
to a lower limit in the metallicity being slightly sub-solar 
(Z ~ 0.8 Zq). This is well within estimates for individual 
young star clusters in the Antennae having metallicities of 
roughly solar (0.5 < Z/Z-. < 1.5. sec Figure 11 in Bastian 



et al.|[2009 l. Metal abundances determined in the hot ISM 
phase are found to vary widely from sub-solar (~ O.2Z0) 
to highly super-solar ~ 20Zq (Fabbiano et al. 2004 Baldi 
et al.||2006a l. So a lower metallicity limit of Z ^ 0.8 Zq is 
not ruled out by the observations. 

On the other hand, we can try and estimate the ob- 
served gas-to-dust ratio. These estimates depend critically 
on the adopted CO-to-H 2 conversion factor. If we identify 
the total gas mass in the disks as M t 



A/4; sk + Mi 



H 2 



( M H 2 



M- 



disk 



since there is no H2 detected in the tidal 

_670 for 

4.78 M pc _/ (K kms -1 ) -1 (|Gao et al .[ 
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quently adopted in starbursts (Downes & Solomon 



and 



2001)1, 

)~ as fre- 



and consistent with the latest analysis (Bolatto et al 



1998) 



2013) 



and values obtained in the nucleus of NGC4038 and the 
overlap regions in the Antennae ( Zhu et aL"||2003 |. Here we 



have used Mi 



H 2 



Qco ■ ico with Leo taken from 



(2001), MSr from 



Hibbard et al. 



Gao et al. 



(2001), and all quantities 



scaled to the same distance. While the uncertainty is large 
depending on our choice for aco, these values seem to favour 
a higher gas-to-dust ratio than the one we have chosen to 
match the observed FIR SED. 

What else could influence the simulated SED shape? 
Changes in the assumed distance used to match the obser- 
vations would result only in changes of the normalization in 
the simulated fluxes, e.g. a distance change of ~ 10% changes 
the fluxes by ~ 20%. They would, however, have no influ- 
ence on the SED shapes. Higher dust masses due to higher 



initial gas fractions or due to cooling and subsequent accre- 
tion from a hot galactic halo (M oster et al.|20Tlj , however, 
might lead to redder simulated SEDs by increasing the total 
dust mass at any given time throughout the merger process, 
without invoking low gas-to-dust ratios. Accordingly higher 
SFRs might counteract this trend, most likely, however, to a 
much smaller extent (see e.g. Hay ward et al.|2011 1. We plan 
to test this hypothesis in future hydrodynamic simulations 
including the effects of metal evolution and the presence of 
galactic hot haloes. 

The specific form of the stellar initial mass function 
(IMF) also affects the simulated SEDs. The star formation 
and feedback model we use ( Springel & Hernquist 2003 1 in- 
trinsically assumes a Salpeter IMF, e.g. for computing the 
supernova rates and hence the feedback energy and gas re- 
turn fractions to the surrounding ISM, producing only about 
half of the luminous stars with respect to a |Chabrier| ( |2003"| 
IMF. For consistency, we also use a Salpeter IMF to com- 
pute the stellar emission in the RT calculations. Choosing 
a Chabrier IMF to test the impact of the IMF in the RT 
calculations, we found that a lower mass-to-light ratio (by 
~ 0.1 — 0.2 dex) in this case leads to an increase in the 
FIR emission by a similar amount (20 — 50%). Further- 
more, the higher relative fraction of massive young stars 
leads to higher dust temperatures, with the emission peaking 
at shorter wavelengths and shifting the SEDs towards the 
blue. As a result, an even lower gas-to-dust ratio would be 
required to match the observed SEDs in Figure [8] Adopting 
a Chabrier IMF in the star formation and feedback model 
would probably yield a more vigorous stellar feedback due 
to the higher fraction in massive stars. The consequences of 
these changes, however, are less obvious, requiring further 
self-consistent simulations. 

The radiative transfer calculation we use here assumes 
that the dust grains are in thermal equilibrium with the 
local radiation field. However, stochastic heating of very 
small dust grains by absorption of starlight results in high- 
temperature transients, during which the absorbed energy is 
re-radiated mainly at mid-infrared wavelengths. It has been 
shown that this mechanism may have a significant effect on 
the emission at 70 /im, depending on the local stellar radia- 
tion field ( |Silva et al.||1998| |Li fe Draine|2001| |. At the FIR 
wavelengths, which are the focus of this article, the emission 
from large grains at the equilibrium temperature, however, 
is expected to be dominant. 

Another factor affecting the results from the radiative 
transfer simulations is the treatment of small-scale structure 
in the ISM. While we assume a smooth density distribution 
at the scale of individual cells (~ 40 pc in the central parts of 
the galaxies), the ISM definitely shows structures on much 
smaller scales. We do not include any sub-scale modelling 
of photodissociation regions (PDRs) around young stars or 
clumpy molecular clouds. In the simulations discussed in this 
article the effect of PDRs is likely to be relatively small, be- 
cause a large fraction of the total luminosity is from stars 
older than 10 Myr. On the other hand, the non-uniform dust 
density distribution is known to have a potentially signifi- 
cant effect on the dust temperature distribution and infra- 



red emission (e.g. Witt & Gordon 1996 Doty et al. 2005 



Schunck et al. 20071. A proper treatment of the clumpy 



cloud structure using a sub-grid model is to be considered 
in future studies. 
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In this paper we have shown as a proof of principle 
that variations in parameters governing the stellar feedback 
model result in unambiguous changes in the observable FIR 
properties. Without attempting a fine-tuned match, in our 
adopted star formation and feedback model a rather inef- 
ficient thermal stellar feedback formulation (0.01 < (?EoS <■ 
0.1) seems to be clearly preferred both by simulated FIR 
maps and SEDs, similar to the recently advocated range by 



Hopkins fe Quataert| ( |2010| ). Recently, more realistic numeri- 
cal models also include metal evolution, kinetic and thermal 



stellar feedback (e.g. Hopkins et al. 2011 1, feedback from 
active galactic nuclei (e.g. Choi et al. 20121, effects from 
magnetic fields (Kotarba et al. 20101 and accretion from 



hot galactic gas haloes (Sinha & Holley-Bockelmann 2009 
Moster et aT1|2011| ). Tests with other codes, e.g. AMR grid 



codes, might also prove useful to gauge numerical artefacts 
inherent to different methods using the same initial condi- 
tions (see e.g. Agertz et al.|200~ for a comparison between 
SPH- and grid-based hydrodynamical methods). The next 
step is to compare further mock observations of the simula- 
tion results to multi-wavelength observations of the Anten- 
nae. In this context, constraints from detailed X-ray maps 



are particularly important (Fabbiano et al. 2001 2004 Baldi 
et al.|p006alb| ). In doing so, using the Antennae system as 
a "standard ruler" , we can further gauge the employed star 
formation algorithms and simultaneously possibly gain fur- 
ther insights into the physics of merger induced star forma- 
tion. 
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